function [f,x,xp,y,yp,symparams,Phi,h1,nomSDF_cup] = DSGEmodel(unitFree)
% This function reports the equations for the New Keynesian model with Calvo prices
% in f, and it reports the state vector (x and xp) and the control vector (y and yp)

% IMPORTANT: lagged controls leaded one period forward must be written as
% the current controls. That is c_cu and not c_ba1p

%Define parameters
syms BETTA B CHI CHI0 THETA DELTA ALFA PHI PHIzero ZI KAPAw ETA RHOR PHIpai PHIpai_1 ...
     PHIy PHIy_1 PHIc PHIc_1 PHIl NU U0 U0d...
     RHOz  RHOd   RHOn  RHOp RHOa ...
     Kss OUTPUTss PAIss MUZss Css AA lss Wss;
Rss = sym('Rss');
symparams=[BETTA,B,CHI,CHI0,THETA,DELTA,ALFA,PHI,PHIzero,ZI,KAPAw,ETA,...
     RHOR,PHIpai,PHIpai_1,PHIy,PHIy_1,PHIc,PHIc_1,PHIl,NU,U0,U0d,...
     RHOz ,RHOd   RHOn  ,RHOp,RHOa, ...
     Kss OUTPUTss PAIss MUZss Css AA lss,Rss,Wss];

%Define variables 
syms c_cu   r_cu    pai_cu   w_cu  l_cu  output_cu  mc_cu  evf_cu  vf_cu  p1_cu  p1Real_cu  dc_cu  ...
     c_cup  r_cup   pai_cup  w_cup l_cup output_cup mc_cup evf_cup vf_cup p1_cup p1Real_cup dc_cup ...
     c_ba1  muz_cu   d_cu  n_cu  paistar_cu  a_cu...
     c_ba1p muz_cup  d_cup n_cup paistar_cup a_cup;
 
%% The equations in f
% EQ1: The expression for the value function (minus)
%vf_cup = - (U0 + d_cup*(1/(1-CHI)*((c_cup-B*c_cu/muz_cup)/Css^CHI0)^(1-CHI)+U0d + n_cup*PHIzero/(1-1/PHI)*(1-l_cup)^(1-1/PHI)) ) + BETTA*AA*evf_cup^(1/(1-ALFA));    
if unitFree == 1
    f1 =  -1 + (-(U0 + d_cu*(1/(1-CHI)*((c_cu -B*c_ba1/muz_cu)/Css^CHI0)^(1-CHI)+U0d + n_cu*PHIzero/(1-1/PHI)*(1-l_cu )^(1-1/PHI)) ) ...
                + BETTA*AA*evf_cu ^(1/(1-ALFA)))/vf_cu;
else
    f1 =  -vf_cu - (U0 + d_cu*(1/(1-CHI)*((c_cu -B*c_ba1/muz_cu)/Css^CHI0)^(1-CHI)+U0d + n_cu*PHIzero/(1-1/PHI)*(1-l_cu )^(1-1/PHI)) )  + BETTA*AA*evf_cu ^(1/(1-ALFA));
end

% EQ 2: Household's First-order condition for labour
if unitFree == 1
    f2 = -1 + (KAPAw*Wss+(1-KAPAw)*(n_cu*PHIzero*(1-l_cu)^(-1/PHI)/((c_cu -B*c_ba1/muz_cu)/Css^CHI0)^(-CHI)*Css^CHI0))/w_cu;
else
    f2 = -w_cu + KAPAw*Wss+(1-KAPAw)*(n_cu*PHIzero*(1-l_cu)^(-1/PHI)/((c_cu -B*c_ba1/muz_cu)/Css^CHI0)^(-CHI)*Css^CHI0);
end

% EQ 3: Euler-equation for the one-period interest rate
% The real stochastic discount factor
mReal_cup = BETTA*d_cup/d_cu*...
    (AA*evf_cu^(1/(1-ALFA))/(vf_cup*muz_cup^((1-CHI)*(1-CHI0)) ))^ALFA...   
    *(c_cup-B*c_cu/muz_cup)^(-CHI)/((c_cu-B*c_ba1/muz_cu)^(-CHI))*muz_cup^(-CHI*(1-CHI0)-CHI0);  
if unitFree == 1
    f3 = -1 + mReal_cup/pai_cup*r_cu;
else
    f3 = -1 + mReal_cup/pai_cup*r_cu;
end

% EQ 4: Firms FOC for labour
if unitFree == 1
    f4 = -1  + (w_cu/((1-THETA)*a_cu *Kss^THETA*l_cu ^-THETA))/mc_cu;
else
    f4 = -mc_cu  + w_cu/((1-THETA)*a_cu *Kss^THETA*l_cu ^-THETA);
end

% EQ 5: First-order condition for prices
if unitFree == 1
    f5 = ((1-ETA)+ZI*mReal_cup*(pai_cup/PAIss^NU-1)*pai_cup/PAIss^NU*output_cup/output_cu*muz_cup ...
        - ZI*(pai_cu/PAIss^NU-1)*pai_cu/PAIss^NU)/(ETA*mc_cu) + 1;
else
    f5 = (1-ETA)+ZI*mReal_cup*(pai_cup/PAIss^NU-1)*pai_cup/PAIss^NU*output_cup/output_cu*muz_cup ...
        + ETA*mc_cu - ZI*(pai_cu/PAIss^NU-1)*pai_cu/PAIss^NU;
end

% Taylor rule
f6= -r_cu + Rss*exp((PHIpai*(log(pai_cu)-(log(PAIss) + log(paistar_cu)) )...
         + PHIpai_1*(log(pai_cup)-(log(PAIss) + log(paistar_cup))) ...
         + PHIy*log(output_cu/OUTPUTss)...
         + PHIy_1*log(output_cup/OUTPUTss)...
         + PHIc*(c_cu-c_ba1+log(muz_cu)-log(MUZss)) ...
         + PHIc_1*(c_cup-c_ba1p+log(muz_cup)-log(MUZss)) ...
         + PHIl*(log(l_cu/lss))));

    
%EQ7: Agg. ressource constraint
if unitFree == 1
    f7 = - output_cu*(1-ZI/2*(pai_cu /PAIss^NU-1)^2)/(c_cu  + DELTA*Kss) + 1;
else
    f7 = - output_cu*(1-ZI/2*(pai_cu /PAIss^NU-1)^2) + (c_cu  + DELTA*Kss);
end

% EQ 8: The aggregated production function to get output. 
if unitFree == 1
   f8 = -1 + a_cu*Kss^THETA*l_cu ^(1-THETA)/output_cu;   
else
   f8 = -output_cu + a_cu *Kss^THETA*l_cu ^(1-THETA);  
end

% The expected value of the value function 
if unitFree == 1
    f9 = -1 + (vf_cup*muz_cup^((1-CHI)*(1-CHI0))/AA)^(1-ALFA)/evf_cu;
else
    f9 = -evf_cu + (vf_cup*muz_cup^((1-CHI)*(1-CHI0))/AA)^(1-ALFA);
end

% Nominal one-period bond price
if unitFree == 1
    f10 = - 1 + mReal_cup/pai_cup/p1_cu;    
else
    f10 = - p1_cu + mReal_cup/pai_cup;    
end

% Real one-period bond price
if unitFree == 1
    f11 = - 1 + mReal_cup/p1Real_cu;    
else
    f11 = - p1Real_cu + mReal_cup;    
end

% consumption growth
if unitFree == 1
    f12 = -dc_cu + log(c_cu) - log(c_ba1) + log(muz_cu);
else
    f12 = -dc_cu + log(c_cu) - log(c_ba1) + log(muz_cu);
end


%% The links
% For lagged consumption
if unitFree == 1
    f13 = -1 + c_cu/c_ba1p;
else
    f13 = -c_ba1p + c_cu;
end

%% The linear shocks
%For non-stationary technology
f14 = -log(muz_cup/MUZss) + RHOz*log(muz_cu/MUZss);

%The preference shock
f15 = -log(d_cup) + RHOd*log(d_cu);

%For labor shocks
f16 = -log(n_cup) + RHOn*log(n_cu);

%Inflation target
f17 = -log(paistar_cup) + RHOp*log(paistar_cu);

%stationary technology shocks
f18 = -log(a_cup) + RHOa*log(a_cu);

%Creating function f
f =[f1; f2; f3; f4; f5; f6; f7; f8; f9; f10;
    f11;f12;f13;f14;f15;f16;f17;f18];

% Define the variables which needs to be log-transformed
y  = [c_cu    r_cu    pai_cu  w_cu    l_cu    output_cu   mc_cu   evf_cu   vf_cu  p1_cu   p1Real_cu  ];
yp = [c_cup   r_cup   pai_cup w_cup   l_cup   output_cup  mc_cup  evf_cup  vf_cup p1_cup  p1Real_cup ];
x  = [c_ba1  muz_cu   d_cu  n_cu  paistar_cu  a_cu];
xp = [c_ba1p muz_cup  d_cup n_cup paistar_cup a_cup];

% Make f a function of the logarithm of the state and control vector
f = subs(f, [x,y,xp,yp], exp([x,y,xp,yp]));
h1 = [];    %endogenous states (not lagged controls)
Phi = [];
nomSDF_cup  = subs(mReal_cup/pai_cup,[x,y,xp,yp], exp([x,y,xp,yp]));

% Define the variables which needs to be log-transformed
y  = [c_cu    r_cu    pai_cu  w_cu    l_cu    output_cu   mc_cu   evf_cu   vf_cu  p1_cu   p1Real_cu  dc_cu];
yp = [c_cup   r_cup   pai_cup w_cup   l_cup   output_cup  mc_cup  evf_cup  vf_cup p1_cup  p1Real_cup dc_cup];
x  = [c_ba1  muz_cu   d_cu  n_cu  paistar_cu  a_cu];
xp = [c_ba1p muz_cup  d_cup n_cup paistar_cup a_cup];




end